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1  Introduction 

The  ocean  is  inhabited  by  innumerable  individuals  of  many  different  genera  and  species, 
each  having  its  own  developmental/physiological  state  and  each  being  immersed  in  its  own 
environment.  The  organisms  move,  both  because  of  water  flow  and  because  of  their  own 
swimming  or  buoyancy,  and  interact  with  their  environment  by  gathering  resources  which 
they  need  and  by  excreting  waste  products.  The  assimilated  material  can  be  used  for 
maintenance,  growth,  or  reproduction.  Finally,  the  organisms  can  die  either  from  natural 
causes  or  because  of  attacks  by  another  organism. 

Furthermore,  the  processes  just  described  must  generally  be  regarded  as  stochastic.  For 
example,  the  probability  of  a  predator  capturing  a  prey  item  will  depend  on  multiple  factors, 
each  with  its  own  probability:  finding  a  prey  item  in  range,  the  choice  to  attack,  success 
in  the  attack,  competition  against  others.  Such  a  description  suggests  an  “agent-based”  or 
“individual-based”  model  (IBM),  with  each  agent  carrying  information  about  its  position, 
its  species,  its  physiological  state,  etc.  Individuals  can  grow,  reproduce,  and  die.  Certainly, 
we  can  build  small  versions  of  such  models,  but  the  number  of  individuals  is  necessarily 
limited  (compared  to  nearly  20,000  copepods  per  cubic  meter  or  to  phytoplankton  densities 
on  the  order  of  108  per  cubic  meter).  However,  such  experiments  may  indeed  give  insight 
into  the  way  in  which  the  local,  stochastic  interactions  translate  into  terms  representing, 
for  example,  grazing  rates  in  terms  of  average  densities. 

Or  we  could  take  the  alternative  view  of  attempting  to  predict  the  probability  distribu¬ 
tion  for  biomass  in  a  continuous  space,  in  this  case  using  something  like  weight  and  species 
as  our  variables.  The  latter  is,  of  course,  discrete,  yet  different  organisms  can  be  genetically 
or,  more  importantly,  functionally  quite  close  to  others.  If  we  choose  a  species  ordering  such 
that  the  maximum  growth  rate  varies  smoothly,  we  may  expect  that  other  terms  entering 
the  dynamics  such  as  the  losses  by  predation  will  also  fall  on  a  fairly  smooth  curve.  Cer¬ 
tainly  on  any  diagram  such  as  figure  (1)  the  gaps  will  be  so  small  as  to  be  negligible,  and 
viewing  the  ordinate  as  a  continuous  variable  is  not  unreasonable.  We  can  then  consider 
the  ways  that  the  processes  described  above  alter  the  biomass  distribution  in  this  space. 

As  in  the  individual-based  model,  the  number  of  variables  we  would  have  to  consider  is 
still  unmanageably  large.  Furthermore,  for  each  (w,  sp ),  we  need  to  specify  the  sources,  the 
sink,  and  the  transfer  rates,  including  possible  nonlinear  dependence  on  the  local  biomass 
6(tc,-s|x,  t)  and  the  density  at  the  source/sink  b(w,s\x,  t). 

Any  attempt  to  construct  such  a  model  or  the  IBM  model,  will  inevitably  point  out 
how  little  we  know  about  most  of  the  species  inhabiting  the  ocean.  We  retreat  to  dealing 
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Figure  1:  The  range  of  sizes  for  various  species 


with  variables  such  as  P  which  can  be  viewed  as  the  integral  over  some  range  (si—  >  s2) 
of  species  and  over  weight.  We  then  presume  that  transfers  into  and  out  of  the  resulting 
black-boxes  can  be  represented  as  functions  just  of  the  integrated  values  and  attempt  to 
parameterize  those. 


2  Basic  biological  models 

Exponential  growth 


Suppose  each  individual  (or  pair)  reproduces  in  a  time  order  r,  then  the  time  evolution 
of  the  population  would  be  described  as 


b(t.  +  St)  =  [1  +  gSt]b(t.),  g  =  ln2/r 
^ =  gb  =>-  b(t)  =  bQexp(gt). 


(1) 

(2) 


where  b  is  the  population  and  g  is  its  growth  rate.  This  model  predicts  an  unrealistic  ex¬ 
ponential  growth  of  the  population  (see  figure  (2)).  Thus,  some  limiting  factors  should  exist. 


Resource  limitation 


For  example,  one  could  write  a  model  with  the  growth  rate  depending  on  a  resource  R 
which  itself  evolves  with  time: 


ft  = 

—  =  -gRb-X(R-R0) 


(3) 

(4) 
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Figure  2:  USA  population  as  a  function  of  time  showing  exponential  growth  at  initial  stage, 
but  a  decay  of  a  growth  rate  at  later  stages 

Here,  the  population  evolves  with  a  growth  rate  gR  and  a  death  rate  d.  The  resource  de¬ 
creases  due  to  consumption  by  population  and  is  restored  to  R0  with  a  time  scale  1/A.  This 
model  allows  for  a  steady  state  solution  with  R  =  d/g,b  =  —\{d/g  —  R0)/d.  However,  the 
growth  rate  is  unconstrained  with  increase  of  the  resource;  a  more  realistic  growth  model 
would  have  to  limit  it  at  R  — >  oo.  One  possibility  is  to  use  the  form  g^j:- 

Predator-pray 

Assume  that  zooplankton  (Z)  grows  by  consuming  phytoplankton  (P)  and  phytoplank¬ 
ton  grows  due  to  source  of  nutrients  (N)  assumed  to  be  fixed;  both  have  natural  death  rates 
dp,  dz 

gPZ  -  pNP  -  dPP  (5) 

agPZ  —  dzZ  (6) 

This  model  (the  Lotka-Volterra  equations)  gives  a  solution  with  cycles  that  have  a  con¬ 
served  quantity  H  =  agP  —  dzlnP  +  gNZ  +  dplnZ.  The  quadratic  NPZ  model  assumes 
N  =  Np  —  P  —  Z,  where  Np  is  a  total  amount  of  nutrients.  This  model  has  three  steady 
states  if  Np  >  (dz/otg)  +  (dp/g),  with  a  P,Z  non-zero  point  being  a  stable  attractor. 

NPZD 

The  NPZ  model  assumes  all  dead  organisms  or  excreted  material  is  immediately  rem¬ 
ineralized  to  usable  nutrient.  In  contrast,  the  NPZD  model  assumes  that  dead  organisms 
and  unassimilated  phytoplankton  would  contribute  to  a  detrital  pool  that  eventually  be- 


dP 

~dt 

dZ_ 

~dt 
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Figure  3:  Schematic  representation  of  NPZD  model  showing  the  fluxes  of  biomass  and  the 
parameters  controling  them 


comes  a  source  of  nutrients.  A  schematic  representation  of  this  model  is  shown  in  figure  (3), 
and  a  typical  version  of  the  corresponding  equations  (including  sinking  of  detritus)  would 
be: 

W,p  ~  '■VT^F-!'pTfeZ-‘ipF  +  v'(KVF)  <7) 

JjiZ  =  09 p  +  kZZ  ~dzZ  +  V'  (8) 

prD  +  -j-wsD  =  (l-a)g  P  Z  +  dPP  +  dzZ  -  rD  +  V  •  (acV D)  (9) 

Dt  oz  F  +  kz 

WtN  =  rD-^p+v'(*™l  (“) 


where  ^  +  u  •  V  denotes  the  material  derivative  of  the  advected  scalar. 


In  general  case  the  system  equations  have  the  form 


-bi  +  V  •  (u ibi)  =  Bi(x,t |b)  +  V  •  ((k  +  Ki)Vbi) 


Here,  u*  and  Kj  represent  the  motion  relative  to  the  fluid.  The  function  Hj(x,  t|b)  represents 
interactions  between  different  parts  of  the  system.  Some  times  a  condition  ^  Hj(x,  i|b)  =  0 
is  used,  which  assumes  that  the  total  biomass  is  conserved;  this  reduces  the  order  of  the 
system. 

Such  biological  models  require  a  large  number  of  parameters.  For  example,  in  the 
NPZD  model  with  4  variables,  the  parameters  are:  fi.  kp,  kz,  g,dp,  dz,  a,  r  in  addition  to  the 
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biological  movement  term  ws  and  physical  parameters  such  as  n;  together,  they  determine 
the  fundamental  behavior  of  the  system.  Thus,  in  addition  to  the  mathematical  question 
of  the  behavior  of  the  equations  for  different  parameter  values,  modellers  need  to  decide 
on  numbers  which  seem  to  give  “realistic”  solutions.  Since  the  variables  are  not  generally 
directly  observed,  comparisons  with  data  tend  to  be  pretty  qualitative  in  an  case. 

Controlled  (i.e.,  laboratory)  experiments  can  be  a  guid  for  parameter  choicese;  however 
the  behavior  of  the  organisms  might  be  affected  by  the  laboratory  constraints  so  that  things 
like  predation  rates  may  be  quite  different.  Another  way  is  to  observe  the  biosystem  in  the 
nature,  but  in  reality  the  system  is  much  more  complicated  than  the  model,  and  this  can 
also  introduce  large  uncertainties. 

3  Stochastic  Dynamics 

An  alternative  method  for  modeling  organism’s  population  dynamics  is  through  the  use 
of  probability  density  functions  for  the  population  size.  The  probability  density  functions 
will  be  written  as  V(n,t)  which  represents  the  probability  of  having  a  population  of  size 
n  at  time  t.  In  case  of  multiple  interacting  populations,  the  evolution  of  the  population 
probability  is  based  on  the  equation 

V(z,  t  +  St)  =  Vt{Sz |z  —  Sz)V(z  —  Sz,  t )  (13) 

where  z  is  a  vector  representing  the  population  of  the  various  species  considered  in  the 
model  and  Vt  is  the  transition  probability.  This  equation  is  essentially  an  application  of 
the  law  of  total  probability.  The  probability  of  a  specific  population  at  a  future  time  step 
is  the  sum  of  all  possible  previous  populations  times  the  probability  of  transitioning  to 
the  specified  population  during  the  given  time  step.  For  one  of  the  species,  the  transition 
probability  can  be  modeled  as 

{ndn5t  Sn  =  —  1 

1  ~  ngnSt  —  ndnSt  Sn  =  0 
ngnSt  Sn  =  1 

where  gn  and  dn  are  the  growth  and  death  rates  per  capita  depending  on  the  population 
size  n. 

To  make  this  more  tractable,  let  us  consider  a  case  where  only  one  population  is  modeled. 
For  simplification  we  will  change  the  notation  from  V(n ,  t)  to  be  Vn(t).  Equation  13  becomes 

Vn{t  +  St)  =  Vn-i{t)(n  -  l)gn-\St  +  Vn(l  -  ngnSt  -  ndnSt)  +  Vn+i (n  +  l)dn+1St 
^nly  +8t>  'Pn  =  'Pn-iitf^n  —  l)g„.-i  —  Vnn(gn  +  dn)  +  Vn+\(n  +  l)dn+\ 

=  Vn-i(n-l)gn-i-Vnnl{gn  + dn) +  Vn+i{n  +  l)dn+i 

O') 

From  this  the  first  few  equations  can  be  calculated  to  give  a  foundation  for  insight 

&Vo  =  Vxdx 

=  -Pi  (91  +di)  +  2V2d2  . 

iV 2  =  Vigi-2r2(g2  +  d2)  +  3V3d3  1  j 

^ \V3  =  2T)2<?2  —  3V3(g3  +  d3)  +  AV3d4 
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Now  we  will  look  at  the  expectation  of  certain  properties  and  will  define  the  functional  form 


^  fm  ^ ^  '  fmVnV 
n 


Clearly  the  sum  of  the  probabilities  at  a  given  time  is  constant  so 


—  <  1  >=  0. 
dt 


A  more  relevant  calculation  is  the  understanding  of  the  time  rate  of  change  of  the  expected 
population  change  which  can  be  shown  to  be 


d 

—  <n  >=<  n(gn 


dn)  >  • 


If  however,  the  death  and  growth  rates  are  independent  of  the  population  size  than  the 
constant  can  be  pulled  outside  of  the  expectation  calculation,  so 


dt<n>=  (9" 


dn)  <n> 


which  has  a  similar  form  to  the  deterministic  population  model.  Additionally  it  can  be 
shown  that  given 


d_ 

dt 


<  n2  >=  2  <  n2(gn 


dn)  >  +  <  n(gn  +  dn)  > 


and  a2  =<  n2  >  —  <  n  >2,  it  is  possible  to  calculate  the  long  term  trend  of  the  ratio 
between  the  standard  deviation  of  the  population  (<t)  and  the  expected  population  size 
with  growth  and  death  rates  independent  of  population  size 


lim  - 

t-yo o  <n> 


N 


o"o  +  <n>  o 

<  n  >o 


g+d 

g-d 


where  do  and  <  n  > o  are  the  initial  standard  deviations  and  the  expected  population  size. 
If  we  compare  a  set  of  realizations,  we  would  find  the  range  of  variation  increases  with  time 
just  as  the  mean  size  is  increasing. 

For  populations  which  can  double  in  a  day,  the  exponential  model  cannot  apply  for  long. 
We  can  examine  the  logistic  model  with  gn  =  go(l  —  n/no),  in  which  case  the  distribution 
can  reach  a  nearly  state.  It  cannot  be  completely  steady,  since  a  constant  value  of  V\  implies 
the  extinction  probability,  Vo,  is  continuing  to  grow. 


3.1  Fokker-Planck  Expansion 

when  the  number  of  individuals  is  very  large,  V(n,  t )  can  be  treated  as  a  smooth  function  of 
n,  regarded  as  a  real  number.  If  we  start  again  with  the  Equation  13,  the  following  Taylor 
expansion  can  be  performed  to  approximate  the  dependence  on  the  vector  of  population 
sizes 

V(z,t)  +  5t^V(z,t)  =  J]^,T(5z|z,t)'P(z,t)  -  Y,SziVT{dz\z,t)V{z,t) 

+  \d^drjYJbzi&zjVT{dz\z,i)V{z,t) 
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Figure  4:  Plot  of  the  U  and  K  for  a  linearly  deacreasing  growth  rate  and  a  constant  death 
rate  and  the  resulting  population  density  distribution. 


where  V(z.  t )  is  still  the  population  probability  function  dependent  on  the  coordinate  z  and 
time,  t ,  and  all  Ts  represent  a  differential  step  in  the  provided  coordinate.  This  Taylor 
expansion  can  be  simplified  to 


d 2 

dzidz-i KtjV 


(16) 


where  the  substitutions 

Ui(z,t)  =  Y  Kij(z,t)  =  Y  VT(Sz\z,t) 

have  been  made.  For  the  case  like  14  considering  birth  and  death  rates  in  one  species, 
U  =  n(gn  —  dn)  and  K  =  \n(gn  +  dn).  An  example  growth  rate  that  can  be  used  is 
9n  =  5o(l  —  n/no)  which  corresponds  to  a  linearly  decreasing  growth  rate  for  populations 
under  a  maximum  population,  uq.  With  a  death  rate  at  a  fixed  constant  the  function 
of  U  and  K  can  be  determined  and  are  presented  in  Figure  4.  The  resulting  population 
probability  distribution  function  is  a  narrow  Gaussian  centered  at  gn  =  dn. 
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